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Abstract. Many complex systems generate multifractal time series which are long- 
range cross-correlated. Numerous methods have been proposed to characterize the 
multifractal nature of these long-range cross correlations. However, several important 
issues about these methods are not well understood and most methods consider only 
one moment order. We study the joint multifractal analysis based on partition 
function with two moment orders, which was initially invented to investigate fluid 
fields, and derive analytically several important properties. We apply the method 
numerically to binomial measures with multifractal cross correlations and bivariate 
fractional Brownian motions without multifractal cross correlations. For binomial 
multifractal measures, the explicit expressions of mass function, singularity strength 
and multifractal spectrum of the cross correlations are derived, which agree excellently 
with the numerical results. We also apply the method to stock market indexes and 
unveil intriguing multifractality in the cross correlations of index volatilities. 
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1. Introduction 

Measurements of a complex evolving system from different view angles provide us many 
time series that are usually long-range cross-correlated and exhibit multifractal nature. 
In turbulent flows, there are velocity held, temperature held and concentration held 
embedded in the same spatial domain. One can measure these quantities at hxed 
locations to obtain time series, which are mutually correlated DEI In hnancial markets, 
there are also many pairs that are cross-correlated, such as market index volatilities, 
price returns of diherent markets, price returns of diherent equities, diherent quantities 
of a same equity [Mn]. Moreover, examples come from very diverse helds, including 
agronomy mm, seismic data ra. meteorology HMS], medical science [la m, 
geophysics [12], transportation pOl - E^] . to list a few. 

To extract the joint multifractality between a pair of multifractal time series, a 
variety of methods have been developed, such as the MF-X-PF method that performs 
joint multifractal analysis [U |2l I23H2S] based on the partition function approach [27] , the 
MF-X-DFA method that conducts multifractal detrended cross-correlation analysis [28] 
based on the detrended huctuation analysis [221 EU], multifractal detrended huctuation 
analysis [3TH33], and the detrended cross-correlation analysis [U EH-HO] , the MF-X-DMA 
method [H] that carries out the multifractal detrended cross-correlation analysis based 
on the detrending moving-average analysis [I2ti52] and multifractal detrending moving- 
average analysis [53l IM]. the multifractal height cross-correlation analysis (MF-HXA) 
method [22], the multiscale multifractal detrended cross-correlation analysis (MM- 
DCCA) [22], and the MF-DPXA method [27] that generalizes the detrended partial 
cross-correlation analysis [2Sl 122] in which the partial correlation is considered. Properly 
designed statistical tests can be used to quantify these cross correlations [601162] . 

The joint multifractal analysis is a classic method and has been applied to study 
the joint multifractal nature between different pairs of time series recorded in natural 
and social sciences [smuEiiiiiiisiEiEn]. Due to its elegant geometric nature, many 
important properties can be derived, which is however very difficult in the frameworks 
of other methods mentioned above. For instance, although there is numerical evidence 
and analytical results for the relationship between the cross-multifractal spectrum 
fxy{o:x,ay) and the multifractal spectra /x(«x) and fy{(yy) of individual time series 
[26l|28|l35] miins], the problem is not solved. Moreover, the original MF-X-PF method 
is important because it handles moments with two different orders, while recent methods 
for multifractal cross-correlation analysis focus only on one order. 

In this work, we recover the uni-order MF-X-PF method [22] and propose a direct 
determination approach for the multifractal spectrum using the idea from the bi-order 
MF-X-PF framework [2]. Based on this framework, we are able to derive important 
geometric properties of the uni-order MF-X-PF method. We perform numerical 
simulations using different mathematical models and explain the results of multifractal 
binomial measures analytically. Finally, we apply the bi-order MF-X-PF method to 
stock market indices. 
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2. Joint multifractal analysis based on partition function approach 

In this section, we first present the joint multifractal analysis based on partition function 
approach with two moment orders [2], abbreviated MF-X-PF(p, g), and then derive the 
uni-order method MF-X-PF(g) that was independently proposed recently [2E]- Although 
the joint partition function XxyiQ, s) of the uni-order method can be directly recovered 
from the joint partition function Xxy{p, Q, s) of the bi-order method by posing p = q, we 
will show that the nexus between the multifractal properties of the two methods is not 
obvious, which is caused by the application of the steepest descent approach. 

2.1. MF-X-PF{p, q) 

Based on the box-counting idea, the geometric support is partitioned into boxes of size 
s. We consider two integrated measures mx{s,t) and my{s,t) in the t-th box. The local 
singularity strengths ax and ay are defined according to the following relationships: 

(la) 

and 

my{s,t) ^ s°‘^. (lb) 

Let Ns{ax, ay) denote the number of boxes of size s needed to cover the set of points in 
which the singularity strengths are around ax and ay with bands dax and day. Hence, 
the fractal dimension of the set is determined according to [6l] 

Ar,(a^,ay) ~ (2) 

in which fxyictxiCty) is the joint distribution of the two singularity strengths |2] or the 
joint multifractal spectrum. 

We consider the joint partition function 

Xxy{p, = [mx{s, t)Y^^ [my{s, (3) 

t 

This dehnition is slightly different from that in Ref. [2], in which the orders are p and 
q rather than p/2 and q/2. In this setting, we recover the traditional partition function 
when mx = rOy and p = q [27]. The joint mass exponent function Txy{p,q) can be 
obtained from the following relation 

Xxy{p,q,s) ^ (4) 

In practice, for a given pair (p, g), we compute Xxy{p,q,s) for a various of box sizes s 
and perform linear regression of \nXxy{p,q,s) against Ins in a proper scaling range to 
obtain Txy{p,q). 

We insert the two relations in Eq. ([^ into the joint partition function, rewrite the 
sum into a double integral over ax and ay, and then apply the steepest descent approach 
to estimate the integral at small s values, which leads to 


'Txyip, q) = P^xl^. + gay /2 - fxy{ax, ay), 


(5) 
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where 

dfxy{ax, OLy)ldax = p/2, (6a) 

and 

dfxy{oix, ay)/day = q/2. (6b) 

Taking partial derivative of Eq. ([^ over p, we have 

dTxy{p,q)/dp = ax/2. (7) 

Similar derivation can be done over q and one can obtain the donble Legendre transforms 


ax = 2dTxy{p,q)/dp, 
ay = 2dTxy{p,q)/dq, 

fxy{ax, ay) = pax{p, q)/2 + qay{p, q)/2- Txy{p, q). 


(8a) 

(8b) 

(8c) 


Therefore, after obtaining Txy{p,q), we can nnmerically determine ax using Eq. (8a), 
ay using Eq. (8b), and fxy{otx,ay) using Eq. (8c). 

From the canonical perspective, one can obtain the fxy{c(x,cty) function directly 
[21 ESI EH]- Dehning the canonical measures as follows 


h'xyiPi S, t) 


[mx{s,t)]P/‘^[my{s,t)]‘i/^ 

^Jm,,(s,f)]p/2[mj/(s,f)]'?/2’ 


(9) 


the two singularity strengths ax{p,q) and ay{p,q) and the joint multifractal spectrum 
fxy{p,q) can be computed by linear regressions in log-log scales using the following 
equations: 

Y.t Pxy{p, q, S, t) Inmxis, t) 


ax{p,q) = lim 


s^o In s 


cty{p,q) = lim 


Et/^xy(p,g,s,f) lnmj,(s,f) 


s-s-o In s 


fxy(,P,q) = lirn 


s^O 


Et Pxyjp, q, S, t) In [Pxyjp, q, S, t)] 
Ins 


(10a) 

(lOb) 

(10c) 


The joint mass exponent function can be obtained by using Eq. ([^. 


2.2. MF-X-PF{q) 

The multifractal cross-correlation analysis based on statistical moments (MFSMXA) 
proposed in Ref. [2H] is actually a special case of MF-X-PF(p, q) when p = q. We call it 
MF-X-PF(q') here for consistency. In this case, we have 

Xxy{q, s) = X] [mx{s, t)my{s, t)^^^ ~ 
t 

in which Txy{q^q) = Txy{q)> Applying the method of steepest descent, Eq. ([^ becomes 

'^xy{_Q') “ 1 “ ^y)l^ fx-yi^^x") 
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where 


d fxyij^xi d fxyic^xi ^y) / d(y.y Q'/2. 


Taking derivative of Eq. (12) over q and using Eq. (13), we have 
dr^yiq) 




q da^ 


a,, 


dq 

Defining that 


= _ + jL--I-^ + - 

2 2 dq 2 ^ 


q doiy h) fxy do^x Idfxy d(y.y oix ^ 

2 dq dax dq day dq 2 


ocxy = [axiq) + ay{q)]/2, 


Eq. (14) and Eq. (12) can be rewritten as follows 

axy = dTxy{q)/dq, 
fxy{oixy{q)) — qo'xy{q) Txy{q), 


(13) 

(14) 

(15) 

(16a) 

(16b) 


where fxyio^xy) — fxy{c(x, cuy). We notice that Eq. (16) has the same form of the 
Legendre transform 

Because ax = dTx{q)/dq and ay = dTy{q)/dq [2Zj, it is easy to verify that the 
following relationship 

'^xyiq) = [Tx{q) + 'ry{q)]/2 + C (17) 

satisfies Eq. (16a), where C is a constant. According to Eq. 10. we have x(0, s) ~ 
sT 2 .h(o) ^ s~^°, which is used to measure the fractal dimension of the geometric support. 
It follows that 

Txy{0) = -Do = -1. (18) 


Combining Eqs. (17) and (18) and using ra;(0) = ry(0) = —1, we have (7 = 0 and thus 

^xy{q) = [rx{q) + ry{q)]/2. (19) 

Inserting Eq. (15) and Eq. (19) into Eq. ( |16b ), we obtain that 

fxy{q) = [fx{q) + /y(g)]/2. (20) 


We note that Eq. (19) and Eq. (20) still hold when Dq ^ 1. In this case, we use 
7 j:2;(0) = T3;(0) = Ty{0) = —Dq to couduct the derivation. These relations were observed 
numerically using the MF-X-DEA method [28], the MF-X-DMA method [H] and the 
MF-X-PF(g) method 


As shown in Eq. (11), the problem is to handle a measure [mx{s,t)^y{s,t)Y^‘^. 
From the canonical perspective, we can obtain the fxyio^xy) function directly 
We can define the canonical measures 


h'xy (*?) S, t) 


[mx{s,t)my{s,t)Y/‘^ 

Y.f(mx{s,t)my{s,t)]<il^' 


( 21 ) 
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The two singularity strengths cxx{p) and cxx{p) and the joint multifractal spectrum 
fxy{p,(l) can be computed by linear regressions in log-log scales using the following 
equations: 


^ \ _ r™ Si Pxyiq, S, t) In [mx{s, t)my{s, _ ax{q) + ay{q) 


where Eq. (10a) and Eq. (10b) are used in the second equality, and 

St Pxy{q,s,t) In [p.xy{q,s,t)] 


fxyi^xyiq)) = lim 


s^O 


In s 


The joint mass exponent function can be obtained by using Eq. (16b). 


(22a) 


(22b) 


3. Joint multifractal analysis of binomial measures 

3.1. Numerical analysis applying MF-X-PF{p,q) 

We perform joint multifractal analysis numerically of two binomial measures [6^. We 
use Px = 0.3 and Py = 0.4 and generate two binomial measures of length 2^°. Figure [^a) 
shows on log-log scales the dependence of Xxy{p,q,s) against box size s for different q 
with fixed p = 2. It is obvious that the curves for different q exhibit excellent power law 
relationships. The power-law exponents obtained by linear regressions of \nxxy{p, q, s) 
against Ins are estimates of the mass exponents Txy{p, q), whose contour plot is shown in 
Fig. [^e). We hud that Txy{p, q) increases with p and q. Adopting the double Legendre 
transform in Eq. (|^, we obtain numerically the singularity functions ax{p^ q) and ay{p, q) 
and the multifractal spectrum fxy{p,q), whose contour plots are illustrated in Fig. St- 
h) respectively. We find that ax{p,q) and ay{p,q) are decreasing functions of p and q, 
while fxy{p,q) has a saddle shape. An intriguing feature is that the contour lines are 
parallel to each other for ax{p, q), ay{p, q) and fxy{p, q)- Figure [^i) plots the singularity 
spectrum fxyictx, ay), which is not a surface but a curve. 

We also calculate the multifractal functions using the direct determination approach 
presented in Eq. (10) for comparison. In Fig. [^b) to Fig. [^d), we illustrate respectively 
the linear dependence of '^^yixy{‘^,q,s,t)\Yi[mx{s,t)], '^^p,xy{2,q, s,t)\n[my{s,t)] and 
J2tPxy{‘^,q,s,t)\n[fixy{,‘^,q,s,t)] against Ins for different q with fixed p = 2. The 
singularity strength functions ax{p,q) and ay{p,q) and the multifractal spectrum 
fxy{p,q) are computed from the slopes of the lines in these three plots. The 
corresponding contour plots are presented in Fig. [^j) to Fig. [^1), which are the same 
as those in Fig. [^f) to Fig. [^h). The numerical results presented in to Fig. [^can be 
derived analytically. 


3.2. Analytical results for MF-X-PF{p, q) 

Let us start with two multifractal binomial measures of length 2^. Consider two 
integrated measures mx{s,t) and my{s,t) in boxes of size s = 2h There are n types 
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Figure 1. Joint multifractal analysis of two binomial measures with = 0.3 and 
Py = 0.4 based on the bi-order MF-X-PF(p, q) method, (a) Power-law dependence of 
Xxy{p,q,s) on box size s for different q with fixed p = 2. (b) Linear dependence of 
s,t)\Ti[mx{s,t)] against Ins for different q with fixed p = 2. (c) Linear 
dependence of l^xy{2., q, s, t) ln[my(s, t)] against Ins for different q with fixed p = 2. 
(d) Linear dependence of /ia;y(2, g, s, t) ln[p.a;j,(2, g, s, t)] against Ins for different q 
with fixed p = 2. (e-i) Mass exponent function Txy{p,q), singularity functions ax{p,q) 
and ay{p,q), and multifractal spectra fxy{p,q) and fxy{oix,Oiy) obtained from (a), 
(j) Singularity function ax{p,q) obtained from (b). (k) Singularity function ay{p,q) 
obtained from (c). (1) Multifractal spectrum fxy{p,q) obtained from (d). 


of boxes whose integrated measures are different, in which 

T 1 T ® 

n = L — I = L — -—. 

In 2 

For the t-the box, we have 

mx{s, t) = = pI{ 1 - PxT~’', 

my{s, t) = t) = 

where k G {1, It follows that 

^ \n my {s, t) — n — Py) 

\npy- Py) 

Inserting Eq. (25) into Eq. ( 24a| ) , we get 

mx{s,t) = C{s) [my{s,t)f = [my{s,t)f, 


/9 = 


\npx - ln(l -px) 
\npy - ln(l -py) 


(23) 

(24a) 

(24b) 

(25) 

(26) 


where 


(27) 
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and 


7 = /?ln(l -p,) - ln(l -p^). 


(28) 


Note that ft and 7 depend only on p^ and Py. When Px + Py = 1, we have (3 = —1. 
When Px = Py, we have /3 = 1 and C'(s) = 1. When both px and Py are greater than 
0.5 or less than 0.5, that is, {px — 0.5){py — 0.5) > 0, we have ft > 0; Otherwise, when 
{Px — 0.5){py — 0.5) < 0, we have ft < 0. 

Combining Eq. ([^ and Eq. (26), we obtain 


Xxyip, q, s) = [my{s, t)] 


Q 


(29) 


where 


Q = (3p/2 + q/2. 


Because roy is a multifractal measure, we have 


^[m^(s,t)] 


^ rs. s'^y(Q) 


(30) 


(31) 


where Ty(Q) has an analytical expression [27j : 

Ty(Q) = -ln[pj + {l-py)^]/\n2. 

The joint partition function can be rewritten as follows 

Xxy{p, q, s) ~ 


Comparing Eq.(|^ and Eq. (33), we obtain the joint mass exponent function: 
7 7 PX ^ 7^7 PI + (1 - Py)'5] 

^xy{p,q) = + = 


2 In 2 


In 2 


It follows that 


- 


2dTxy{p,q) 7 /S P^^ripy +{l-py)^ln{l-py) 


Py + (1 - Py)^ 


and 


ay 


dp In 2 In 2 

2dTxy{p,q) _ 1 p^lnpy + {l-py)Qln{l-py) 


dq In 2 pQ + (i_p^)Q 

We obtain immediately the relationship between and ay 

7 

^ pay. 

In 2 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 


This relationship explains the observation in Fig. j^i) that fxy{c(x,cty) is a curve along 
this line rather than a surface and the line segment ( [^ is the projection of fxy{c(x, oty) 
onto the {ax, ay) plane. 
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We now derive the main geometric properties of ax{Q) and ay{Q). We find that 


cty{Q) is a monotonically decreasing function of Q, because 


day 

dQ 


1 [Inpj; - ln(l 


In 2 


Vy + 


< 0 . 


(38) 


We can prove that the limits of ay exist when Q —)■ ±oo. We rewrite Eq. (36) as follows 

1 \nPy+[{l-Py)/Py]^\n{l-Py) 


We can obtain that 


In 2 


o^y,mm = hm tty = miu 

’ Q—^oo 


tty,max = lim ay = max 

Q—)■—oo 


1 + [{l-Py)/Pyf 

Inpy ln(l -Py) 

In 2 ’ In 2 
\npy ln{l-py) 


(39) 


ln2 ’ 


In 2 


(40) 


Therefore, the solution of Eq. (36) exists and is unique if and only if ay G «y,r 

The explicit form of the solution is 


Q = In 


\0g2[{l -Py)/Py 


- 1 


ay+ \Og 2 Py 

Further, the width of the singularity spectrum of ay is 

|ln(l -py) - \npy\ 


/In 

Py 

/m 

.i^-Py). 


(41) 


Aay = 


In 2 


(42) 


These results explain the parallel observation of the contour lines in Fig. [^g). When 
Py = 0.5, Aay = 0. In this case, the measure is neither multifractal nor monofractal 
since it is uniformly distributed on the support. 


According to Eq. (37), we have 


da^ 

dQ 


f3 p^{l-Py)^[\npy-\n{l-py)f 


In 2 


P? + i^-Py)^ 


(43) 


which suggests that a^ is a strictly monotonic function of Q. Moreover, it is easy to 
show that 

Inpx ln(l-p3,) 

— llilli S — 

Q^oo 


Oix,mm = Jim ax = min 


ttx.max = lim ax = max 
Q^—oo 


ln2 ’ ln2 
Inpx ln{l-px) 


ln2 ’ 


In 2 


(44) 


Therefore, the solution of Eq. (|35|) exists, which is unique if and only if ax G 


'^x,ramy ^x,max} 


]. Due to the symmetry between the two measures mx and my, the results 
for ax are obvious, provided that we know the geometric properties of ay. 
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We now turn to investigate the geometric properties of the multifractal spectrum 
fxy{c(x, Ciy), whlch has the following form: 

fxy{a^, ay) = pa^/2 + qay/2 - T^y{p, q) 


= ? ( ^ 


+ jJoii 


2 Vln2 ' ^ V ' 2 ^ 




“h 

Q 


VI ^ +(1-Pj/)®] 


21n2 


In 2 


Q (^) “ Py'^ ^ 

In 2 / ^ „ \ Q 


In 


1 + 


i-Pa 

Py 


Q 


Q 


1 + 


^-Py 

Py 


Q 


In 


i-Py 

Py 

Py 

I-Pm 


+ 


1 + 


^-Py 

Py 


Q 


In 


In 2 


1 + 


(45) 


l-Pa 

Py 


Q 


In 2 


1 + 


i-Pa 

Py 


Q 


It is easy to hnd that 

fxyiQ = 0) = 1 and fxyiQ) = fxy{-Q), (46) 

where fxyiQ) — fxyiax, ayiQ). It indicates that fxyiax,ay) is symmetric with respect 
to the line Q = 0, as numerically shown in Fig. [^h). Furthermore, we obtain 

(47) 


lim fxyip,q) = 0- 

Q^±oo 

Taking derivative of fxyiax, ay) with respect to Q, we have 


dfxyiQ) Q 


dQ 


In 2 


In 


Py 


Q-Py) 


Py 


l-Py 


< 3/2 


+ 


1 - 




Q/2' 


Py J 


-2 


(48) 


When Q < 0, dfxyiQ)/dQ > 0 so that fxyiQ) is a monotonically increasing function of 
Q. When Q > 0, dfxyiQ)/dQ < 0 so that fxyiQ) is a monotonically decreasing function 
of Q. Therefore, the maximum of fxyiQ) is 1 and its minimum is 0. These properties 
explain the parallel feature of the contour lines in Fig. [^h). 

We note that the numerical results are in excellent agreement with the analytical 


results for Txyip.q) in Eq. Q, axip,q) in Eq. (35), ayip,q) in Eq. (36), and fxyip,q) 


in Eq. (45). Combining Eq. (41) and Eq. (45), we hnd that fxyiax, ay) is a univariate 
function of ay, or of ax by using Eq. (37). 


3.3. Numerical analysis applying MF-X-PFiq) 

We also apply the MF-X-PF(g) method to the same mathematical example. The results 
are shown in Fig. We hnd that the three theoretical relationships in Eq. (19), Eq. (15), 


and Eq. (20) are nicely verihed. In addition, we observe again that the results from the 


classic partition function approach and the direct determination approach agree with 
each other. We note that this is also the case for other mathematical and empirical 
examples investigated in this work. Thus we will not show the results obtained from 
the direct determination approach in the rest of this paper. 
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Figure 2. Joint multifractal analysis of two binomial measures with = 0.3 and Py = 
0.4 based on the MF-X-PF(g) method, (a) Power-law dependence of Xxy{q,s) on box 
size s for different q. (b) Linear dependence of l^xyiq, s,t)ln[mx{s,t)my{s,t)]^^‘^ 
against Ins. (c) Linear dependence of ^^xy{q, s,t)hi[ij,xy{q, s,t)] against Ins. (d) 
The mass exponent function Txy{q). (e) The singularity strength function a{q). (f) 
The multifractal singularity spectrum fxy{a). 


4. Joint multifractal analysis of bivariate fractional Brownian motions 


We further investigate the MF-X-PF(p, q) algorithm using monofractal measures. If mx 
and my are monofractal, we have ax = oty = 1 and fxy = 1 according to its dehnition in 
Eq. 0 IZHEH]. Together with Eq. ( [^ , we have 


r. 


xy 


(p. 9 ) = p /2 + 9 / 2 - 1 . 


(49) 


These properties are indicators of monofractality. 

The mathematical model used here is bivariate fractional Brownian motions 
(BFBMs). The two components x{t) and y{t) of the BFBM are two univariate fractional 
Brownian motions with Hurst indices Hxx and Hyy, respectively. The basic properties 
of multivariate fractional Brownian motions have been comprehensively studied [69Hn]. 
Extensive numerical experiments of other MF-DCCA algorithms have been conducted 
using bivariate fractional Brownian motions miEzi. The two Hurst indexes Hxx and 
Hyy of the two univariate FBMs and their cross-correlation coefficient p are input 
arguments of the simulation algorithm. By using the simulation procedure described in 
Refs. [701 ITT] , we have generated as an example a realization of BFBM with Hxx = 0.1, 
Hyy = 0.5 and p = 0.5. The length of the BFBM is 2^®. The joint multifractal analysis 
of the BFBM using the MF-X-PF(p, g) algorithm is presented in Fig. 

The corresponding power-law dependence of the joint partition function Xxy{p, Q, s) 
with respect to the box size s for different g’s and hxed p = 2 is shown in Fig. 0a). 
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Figure 3. Joint multifractal analysis of bivariate fractional Brownian motions with 
Hxx = 0-1) Hyy = 0-5 and p = 0.5. (a) Power-law dependence of XxyiP,Qjs) on box 
size s for different q with fixed p = 2. (b) Mass exponent function Txy(j),q) function 
obtained from (a), (c) Errors ATxy{p,q) between the estimated exponent Txy{p,q) 

and the theoretical function p/2 -|- q/2 — 1. (d,e,f) Singularity functions ax{p,q) and 
ay{p,q), and multifractal spectra fxyi,o.x,ciy) obtained from (b). 


The scaling ranges span over two orders of magnitude. The slopes of the lines give 
the estimates of Txy{p,q), where p and q vary from —10 to 10 with a spacing of 0.1. 
The resulting mass exponents Txy{p^ q) are shown in the contour plot of Fig. |^b). We 
observe that Txy{p,q) increases with p and q, the contour curves are parallel lines, and 
the parallel lines are evenly spaced. These features suggest that r^yip^q) is a linear 
function of p and g, which is an indicator of monofractality. 

In order to further show the performance of the MFXPF algorithm, we calculate 
the errors between the estimated exponents Txy{p,q) and the theoretical exponents as 
^Txy{p,q) — Txy {p,q) - \p/2 + q/2 - 1]. Fig. |gc ) shows the dependence of ATxy{p,q) 
with respect to p and q. All the APxyip, q) values are less than 0.15, implying that the 
algorithm gives good estimates. 

By adopting the double Legendre transform in Eq. ([^ numerically, we get the 
singularity strength functions ax{p,q) and ay{p,q) and the multifractal spectrum 
fxy{c(x, cty), whose contour plots are shown in Fig. [^d,e,f). The singularity strength 
functions cxxip.q) and ay{p,q) are close to 1, indicating that the functions ax{p,q) and 
ay{p, q) are independent of the order p and q. Although there is a trend in each function 
ax{p, q) and ay(p, g), the theoretical functions axip, q) = ^ and ay{p, g) = 1 are basically 
satished. Hence, the MF-X-PF algorithm is able to correctly capture the monofractal 
nature of the BFBMs. 

Fig. I^g) plots the singularity spectrum fxy{oix,ay), which is a surface and the 
contour lines are closed curves. It is easy to find that the vast majority of the surface 
































Joint multifractal analysis based on the partition function approach 


13 


is nearly equal to the theoretical function fxy{o.x,Oiy) = 1. We observe that the errors 
ATxy{p, q) is equal to the difference between fxy{p, q) and 1, as shown by the Legendre 
transform. 

We point out that the results using the direct determination approach are exactly 
the same as shown in Fig. We thus summarize that the theoretical analysis is well 
verihed by the numerical results. 

5. Application to stock market indexes 

We now apply the MF-X-PF(p, g) algorithm to investigate the long-range power-law 
cross correlations of the daily volatility time series of the Dow Jones Industrial Average 
(DJIA) and the National Association of Securities Dealers Automated Quotations 
(NASDAQ) index. The daily volatility is dehned as the absolute value of the logarithmic 
difference of daily closing prices: 

R{t) = \\nP{t)-\nP{t-l)\, (50) 

where P{t) is the closing price on day t and has been retrieved for the DJIA and 
NASDAQ indices. The time period of the samples is from 5 February 1971 to 25 
January 2011, containing 10084 data points. The daily return time series of the two 
indexes are shown in Figure SI (New J. Phys. online). 



Figure 4. Joint multifractal analysis of the cross correlations between the daily 
volatility time series of DJIA index and NASDAQ index using the MF-X-PF(p, q) 
approach, (a) Power-law dependence of Xxyip, s) on box size s for different q’s with 
fixed p = 2. (b) Mass exponent function T^y{p,q). (c) Singularity strength function 
otx(p,q)- (d) Singularity strength function ay{p,q). (e) Multifractal function fxy{p,q)- 
(f) Multifractal singularity spectrum fxyictx,cXy). 

Fig-i a) shows on log-log scales the dependence of the joint partition function 
Xxyip ) -s) with respect to the box size s for different g’s and hxed p = 2. We observe nice 
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power-law scaling over about 1.5 orders of magnitude. The contour plot of the exponents 
Txy{p, q) is shown in Fig. |^b), where p and q vary from —10 to 10 with a spacing of 0.1. 
The contour curves are not straight lines and the spacings between neighboring curves 
are not equidistant. Fig. |^c) and Fig. |^d) illustrate respectively the contour plots of 
the singularity strength functions a^ip^q) and ay{p,q), which are obtained numerically 
from Txy{p,q). We observe that the values of the singularity strength range from 0.6 
to 1.2, which are well dispersed. In addition, the singularity strength functions are not 
monotonic with respect to p or q. Fig. |^e) illustrates the multifractal function fxy{p, q) 
obtained from the Legendre transform, whose values range from 0 to 1. The maximum 
fxy{p,q) = 1 is reached at point (p, g) = (0,0). Within the investigated intervals of p 
and q, the small fxy{p, q) values concentrated in the region with large values of p and q. 
In Fig. I^f), we present the singularity spectrum fxy{c(x, Oiy). These empirical hndings 
suggest that the cross correlations between daily volatilities of DJIA and NASDAQ 
possess multifractal nature, which is consistent with previous results using the MF-X- 
DFA, MF-X-DMA and MF-X-PF(g) methods [26l EHl SB 172] . 



CHx CKx 


Figure 5. Comparison of the joint multifractal singularity spectrum fxyio:x, o:y) 
between the daily volatility time series of the DJIA index and the NASDAQ index 
in different time periods with and without market turmoil. 


To reveal whether the joint multifractality between the daily volatilities of the two 
indices remains or changes along time, we perform the MF-X-PF(p, q) analysis in moving 
windows on a decade basis with a step of one year. The results are presented in Figure S2 
(New J. Phys. online). We show six plots in Fig. We find that the joint multifractal 
singularity spectrum fxy{c(x, C(y) changes over time. Moreover, the inclusion or exclusion 
of hnancial turmoils (high volatile periods) has a significant impact on the shape of 
fxy{cix,c(y)- Ih the sample period under investigation, there were two infamous market 
crises, the Black Monday in 1987 and the latest crisis in 2008. During relatively calm 
periods, the fxyic^x, C(y) contour looks roughly like an American football. However, when 








Joint multifractal analysis based on the partition function approach 


15 


one of the crisis is included, the contours are significantly stretched to the southwest. 
In other words, the singularity strengthes and ay have much smaller values during 
turmoil periods. This is actually not surprising because this feature is well-documented 
for ordinary multifractals |27j . 

We repeat the same analysis for two stocks Du Pont (NYSE:DD) and Exxon Mobil 
(NYSE:XOM) over time period from 05-Jan-1970 to Ol-Sep-2015, containing 11522 data 
points. The daily return time series of the two stocks are shown in Figure S3 (New J. 
Phys. online). The results are illustrated in Figure S4 (New J. Phys. online). As 
expected, very similar results are observed. Two more pairs of financial time series are 
investigated and the results are presented in Figure S5 to Figure S8 of the Supplementary 
data (New J. Phys. online). One pair is about crude oil commodities, Arab Light to USA 
and WTI Cushing. The sample period from is 03-Jan-1991 to 18-Dec-2012, containing 
5510 data points. Another pair is about Special Drawing Rights (SDRs) per currency 
unit for the U.K. pound sterling (GBP) and the U.S. dollar (USD) over time period 
from 05-Jan-1994 to Ol-Sep-2015, containing 5452 data points. 

Compared with the results of binomial measures and fractional Brownian motions, 
the multifractal function and the multifractal singularity spectrum exhibit different 
shapes for different data sets studied. For example, in Fig. |^f) for the financial market 
data there is a pronounced asymmetry, and the spectrum exhibits a stretched shape, 
in sharp contrast to Fig. [^f) for the artihcial BFBM data. These features reflect the 
irregular nonlinear traits of hnancial indexes. Roughly, the spectrum contour parallels 
to the diagonal ax = ay (cf. Eq. (37)), which is due to the fact that the DJIA and 
NASDAQ indexes comove along time so that the volatilities fulfill Eq. (26) to certain 
extend. A direct conjecture is that the correlation coefficient p{ax, ay) is greater if the 
correlation coefficient p{Rxit), Ry{t)) is greater. This is validated by Fig. S9 in the 
Supplementary Data (New J. Phys. online). 


6. Conclusions 

We have studied the properties of joint multifractal analysis based on partition function 
with two moment orders, termed MF-X-PF(p, g). The uni-order method MF-X- 
PF(g) has then been derived. The main properties of these methods have been 
obtained analytically. For instance, for the MF-X-PF(g) method, we have obtained the 
relationship between the joint mass exponent function and the individual mass exponent 
functions, Txy{q) = [Tx{q) + 'Ty{q)]/‘2, which was numerically and empirically observed in 
the literature. 

We applied the MF-X-PF(p, g) method to multifractal binomial measures. The 
expressions of mass function, singularity strength and multifractal spectrum of the 
cross correlations have been derived, which agree excellently with the numerical results. 
We further validated the performance of the method by using bivariate fractional 
Brownian motions without multifractal cross correlations. When applied to the daily 
volatility time series of two stock market indexes, intriguing multifractality in the cross 
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correlations is confirmed. The multifractal properties of these examples are found to 
be the same when we use the conventional determination approach and the direct 
determination approach. 

Multifractal cross-correlation analysis has been applied in many helds, especially 
in Econophysics. Although there are numerous methods, most of them consider only 
one moment order. It is natural that bi-order methods such as MF-X-PF(p, g) can be 
developed for other uni-order methods. We expect that such bi-order methods will unveil 
new stylized facts in the analysis of hnancial time series, which can serve to calibrate 
agent-based models [73]. In addition, the joint multifractal nature extracted from two 
long-range cross-correlated time series has potential applications. One possibility is to 
construct a multi-scale cross-correlation measure, analogous to other DCCA coefficients 
[SHI ED IS21 EH ESI- Another possibility is to construct a measure quantifying market 
efficiency [761179] . A related possibility is to quantitatively characterize the degree of 
market unrest other than the volatility measure [80] . 
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